#include "VGFAsympExp.h"
#include <gsl\gsl_sf_gamma.h>
#include <boost\math\special_functions\factorials.hpp>
#include <math.h>

VGFAsympExp::VGFAsympExp( double mu, double lambda, double alpha, double beta, size_t nTerms ) 
	: _mu(mu), _lambda(lambda), _alpha(alpha), _beta(beta), _nTerms(nTerms), _rightCoefficients(nTerms, 0), _leftCoefficients(nTerms, 0)
{
	double gamma = sqrt((alpha*alpha) - (beta*beta));
	_const1 = pow(gamma,2*lambda)/(pow(2,lambda)*pow(alpha,lambda)*gsl_sf_gamma(lambda));

	for(int k = 0; k < _nTerms; k++)
	{
		_rightCoefficients[k] = b(k, _beta);
		_leftCoefficients[k] = b(k, -_beta);
	}
}


VGFAsympExp::~VGFAsympExp(void)
{
}

double VGFAsympExp::a( unsigned int k, double v ) const
{
	using namespace boost::math;
	return (pow(-0.5, (int)k)*gsl_sf_poch(0.5 - v,k)*gsl_sf_poch(0.5 + v, k))/factorial<double>(k);
}

double VGFAsympExp::b( unsigned int k, double beta ) const
{
	double sum = 0.0;
	for(size_t j=0; j<=k; j++)
	{
		double product = 1.0;
		for(size_t i=0; i<j; i++)
		{
			product *= (_lambda-k+i);
		}

		sum += product*(pow((_alpha-beta), -int(j)-1)*a(k-j,_lambda-0.5))*pow(_alpha, int(j)-int(k));
	}
	return sum;
}

double VGFAsympExp::LeftTail( double x ) const
{
	x = _mu - x;
	return evaluate(x, -_beta, _leftCoefficients);
}

double VGFAsympExp::RightTail( double x ) const
{
	x -= _mu;
	return 1-evaluate(x, _beta, _rightCoefficients);
	
}

double VGFAsympExp::evaluate( double x, double beta, const std::vector<double>& coefficients ) const
{
	double sum = 0.0;
	for(int k=0; k < _nTerms; k++)
	{
		sum += (coefficients[k]*pow(x,-k));
	}

	sum *= _const1*exp((-x)*(_alpha-beta))*pow(x,_lambda-1);
	return sum;
}
